% Taking the firm-side calibration as given, this script does the following:
% 1) Given the resitriction on absolute fee level, it finds the highest
% rate of return Rz <= 5% that the providers could offer while remaining on the
% laissez-faire iso-profit curve. 
% 2) It verifies whether or not this contract now prevails in equilibrium.
% 3) What is the impact on firms profits and consumer welfare?


clear
load('recalibrated_data_with_eval.mat')

fileID = fopen('Policy_1_ceiling_2021.txt','wt')


% So, what fraction of the original equilibrium charge do we allow?
fraction = 0.35;

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Policy experiment with a CEILING ON FEES equal to %f of the equilibrium fee \r\n'],[fraction]);


rho = mydata.rho;
delta = mydata.delta;
effhh_size_ = mydata.effhh_size;
avgIncome = mydata.avgIncome;
meanhhs = mydata.meanhhs;
R = mydata.R;
zpen_eval = mydata.zpen_eval;
alpha = mydata.alpha;
survival = mydata.survival;
age_ = mydata.age;

THETAtarget = 0.05;
totalcosttarget = 0.005;
adminratiotarget = 0.50;
markuptarget = 0.2;

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_naif_target = mydata.(suffix);
suffix = strcat('fmax_',num2str(THETAtarget*10000));
fmax_target = mydata.(suffix);


% Set the calibrated values:
gamma1 = 0.5;
c1 = adminratiotarget*totalcosttarget*mean(avgZ_naif_target) / mean(avgZ_naif_target.^(gamma1));
gamma2 = 5.75;
c2 = (1-adminratiotarget)*totalcosttarget*mean(avgZ_naif_target) / (1+THETAtarget-R)^(gamma2);
phi = 0.35;
xi = 0.481482;

% And for later comparisons with laissez-faire:
lf_profit = 15211.219049;
lf_actual_utility = 654.673018;


% Re-calculate the iso-fees (relative to the lf profit level, but taking into account the ceiling)
for THETA = mydata.THETA_grid
    if abs(THETA - THETAtarget) < eps
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = fraction*phi*fmax_target;
    else
        suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
        avgZ_ = mydata.(suffix);
                        
        costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
                        
        fee = 100;
        profits = ones(1,71)*fee - costs;
        ProfitPV = profits(1);
        for t = 2:71
            ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
        end
                    
        while lf_profit > ProfitPV
            fee = fee + 6;
            profits = ones(1,71)*fee - costs;
            ProfitPV = profits(1);
            for t = 2:71
                ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
            end
        end
                        
        suffix = strcat('fmax_',num2str(THETA*10000));
        fmax = mydata.(suffix);
                        
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = min(min(fee-3,fmax),fraction*phi*fmax_target);
    end
end

                    
% Calculate attainable profits for each THETA 
profit_grid = [];
fprintf(fileID, ['**************\r\n']);
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    
    suffix = strcat('isofee_',num2str(THETA*10000));
    fee = workspace.(suffix);
         
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
             
    fprintf(fileID, ['For Rz = %f, the PV of profits (on the iso-profit curve) is %f \r\n'],[THETA,ProfitPV]);
    
    profit_grid = [profit_grid ProfitPV];
end


% Which of the contracts prevails in the market now?
% The best initial guess is to look for a contract with maximum THETA
% that allows the firms to stay on the laissez-faire iso-profit, 
% since the consumers preferences peak at 5%.
max_profit = max(profit_grid);

% This check is built in to account for small differences in profit levels
% due to approximation of iso-fees. It isolates the highest THETA from the
% l-f iso-profit as the candidate solution.
id_max_profit = length(mydata.THETA_grid);
for THETA = flip(mydata.THETA_grid)
    if abs(profit_grid(id_max_profit) - max_profit)/max_profit < 0.05
        break
    end
    
    id_max_profit = id_max_profit - 1;
end

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['The profit-maximising contract now has Rz = %f \r\n'],[mydata.THETA_grid(id_max_profit)]);


% In addition, verify that there is no incentive to price-deviate from
% iso_fee(id_max_profit) by +/- 25%: 
THETA = mydata.THETA_grid(id_max_profit);
suffix = strcat('isofee_',num2str(THETA*10000));
fee_supp_eq = workspace.(suffix);
suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
avgZ_ = mydata.(suffix);

suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
avgZ_TC = mydata.(suffix);
suffix = strcat('avgX',num2str(THETA*10000),'_TC');
avgX_TC = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
avgTC_TC = mydata.(suffix);

if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

% Calculate perceived utility at current eqm candidate
avgTC_fee = avgTC_TC - fee_supp_eq*ones(1,71);    
instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
lifetime_new_eq = instant(1);
for t = 2:71
    lifetime_new_eq = lifetime_new_eq + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

% Verify no profitable price deviation: 
gridlength = 11;
gridjump = 0.05;

ProfitPVvec = zeros(1,gridlength);
expected_profits = zeros(1,gridlength);

for j = 1:gridlength
    % First, profit from an interaction
    fee = min((1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*fee_supp_eq, fraction*phi*fmax_target);
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
    ProfitPVvec(j) = ProfitPV;
    
    % Second, probability of interaction
    avgTC_fee = avgTC_TC - fee*ones(1,71);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    lifetime_fee = instant(1);
    for t = 2:71
        lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    probabilityraw = 0.5 + (lifetime_fee - lifetime_new_eq)/(2*xi);
    probability = max(0,min(1,probabilityraw));
    
    expected_profits(j) = ProfitPV*probability;
end

[max_profit,adjustment_id] = max(expected_profits);
optimal_adjustment = (1 - floor(gridlength/2)*gridjump + (adjustment_id-1)*gridjump);
fee_eq = min(optimal_adjustment*fee_supp_eq, fraction*phi*fmax_target);
fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['For xi = %f and Rz = %f, there is no profitable price-based deviation for fee = %f  \r\n '],[xi,THETA,fee_eq]);
fprintf(fileID, ['(Adjustment (if any) with optimal multiplier of  %f implemented)  \r\n '],[optimal_adjustment]);



% Impact on firms profits:
% Note: Need to multiply max_profit by 1/0.5 to correct for eqm
% probability of interaction.
profitchange = (max_profit*2 - lf_profit)/lf_profit;
fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Relative to laissez faire, the change in firms profits is %f percent\r\n'],[100*profitchange]);


% Impact on consumer welfare:
% Calculate ACTUAL utility from the new equilibrium contract and evaluate
% welfare loss, relative to laissez-faire. 
fprintf(fileID, ['**************\r\n']);
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

actual_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid(id_max_profit)
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_naif');
    avgX_ = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
    avgTC_ = mydata.(suffix);

    fee = fee_eq;
                    
    avgTC_ = avgTC_ - ones(1,71)*fee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    %actual_lifetime_iso_grid = [actual_lifetime_iso_grid lifetime];
    fprintf(fileID, ['For Rz = %f and the associated iso-fee, the ACTUAL lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
end

% What is the CE welfare loss?
utility_multiplied = lifetime;

efficient = lf_actual_utility;
multiplier = 1.00;
while utility_multiplied < efficient 
    instant_multiplied = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
        else
            instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    utility_multiplied = instant_multiplied(1);
    for t = 2:71
        utility_multiplied = utility_multiplied + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
    end
    multiplier = multiplier + 0.0001;
end

fprintf(fileID, ['The CE welfare loss relative to laissez faire is %f percent\r\n'],[100*(multiplier-1.0001)]);








    
    
                    
                    
